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Abstract. Monte Carlo simulations of the two-dimensional XY model are 
performed in a square geometry with various boundary conditions (BC). Using 
conformal mappings we deduce the exponent r?o- (T) of the order parameter correlation 
function and its surface analogue ryy (T) as a function of the temperature in the 
critical (low-temperature) phase of the model. The temperature dependence of both 
exponents is obtained numerically with a good accuracy up to the Kosterlitz-Thouless 
transition temperature. The bulk exponent follows from simulations of correlation 
functions with periodic boundary conditions or order parameter profiles with open 
boundary conditions and with symmetry breaking surface fields. At very low 
temperatures, r\ a (T) is found in a pretty good agreement with the linear temperature- 
dependence of Berezinskii's spin wave approximation. We also show some evidence 
that there are no noticeable logarithmic corrections to the behaviour of the order 
parameter density profile at the Kosterlitz-Thouless (KT) transition temperature, 
while these corrections exist for the correlation function. At the KT transition 
the value »7ct(Tkt) = 1/4 is accurately recovered. The exponent associated to the 
surface correlations is similarly obtained after a slight modification of the boundary 
conditions: the correlation function is computed with free BC, and the profile with 
mixed fixed-free BC. It exhibits a monotonous behaviour with temperature, starting 
linearly according to the spin wave approximation and increasing up to a value 
7?|| (Tjct) — 1/2 at the Kosterlitz-Thouless transition temperature. The thermal 
exponent r\ e (T) is also computed and we give some evidence that it keeps a constant 
value in agreement with the marginality condition of the temperature field below the 
KT transition. 
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1. Introduction 

The two-dimensional classical XY model has a rich variety of applications in condensed 
matter physics . It describes of course magnetic films with planar anisotropy f , but 
also thin-film superfluids or superconductors, or two-dimensional solids. In statistical 
physics, this model was also extensively studied for fundamental reasons, as describing 
for example classical Coulomb gas or fluctuating surfaces and the roughness transition. 
Although no exact solution exists for this model, many of its essential properties are 
known from different approaches. 

The model undergoes a standard temperature-driven paramagnetic to ferromag- 
netic phase transition in d > 2, characterised e.g. by a power-law divergence of the 
correlation length near criticality, 

£~ \t\~ v , (1) 

(see e.g. Ref. Q). In two dimensions, it exhibits a rather different behaviour, 
since long-range order of spins with continuous symmetries (if there is no source of 
anisotropy in spin space) is forbidden according to the Mermin- Wagner theorem ||, . 
Indeed there is no spontaneous magnetisation at finite temperature in two-dimensional 
systems having a continuous symmetry group, for example XY or 0(2) spin models, 
since order would otherwise be destroyed by spin wave excitations [||. On the other 
hand, there is an infinite susceptibility at low temperatures, indicating another type of 
transition ||. Such systems can indeed display a low-temperature phase with quasi- 
long-range order and a defect-mediated transition towards a usual paramagnetic phase 
at high temperatures. This transition is governed by unbinding of topological defects. 

We now specify the case of the 2d AF-model: consider a square lattice with two- 
components spin variables cr w = (cos9 w ,sm9 w ), |<x lu | 2 = 1, located at the sites w of 
a lattice A of linear extent L, and interacting through the usual nearest-neighbour 
ferromagnetic interaction 

k <j w ■ <t w+ p l 



W fl 



= kJ2J2 cos ( 9 ™ - e »+A). ( 2 ) 



W ft 



where K = J/ksT, with J the ferromagnetic coupling strength and fj, = 1,2 labels 
the directions in A (ft is the corresponding unit vector). The angle 9 W measures the 
deviation of tr w from an arbitarty reference direction. 

The low-temperature regime, governed by spin-wave excitations, was investigated 
by Berezinskii || in the harmonic approximation. In the low-temperature limit 
(also called spin- wave phase), localised non- linear excitations (called vortices) are 
indeed bounded in pairs of zero vorticity, and thus do not affect crucially the spin- 
wave description. At a first approximation the effect of vortices can be completely 
neglected and the spin wave approximation, after expanding the cosine, reproduces 
the main features of the low-temperature physics. This harmonic approximation is 
justified, provided that the spin disorientation remains small, i.e. at sufficiently low 



f Id X y-type singularities are extremely hard to observe experimentally in real magnetic layers, 
since they are washed out by crossover effects as soon as a small interlayer interaction takes place in 
the system. 
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temperature: 



H H n 



k B T k B T 



I^EE^-w 2 - (3) 

w ^ 

Within this approximation, the quadratic energy corresponds to a Gaussian 
equilibrium distribution and the two-point correlation function at inverse temperature 
(3 becomes fj, §, §, 

z ? = ( n / im- w« - w) a ]. (4) 



W.f-L 



e -3vnj^(»-.-»«.)'«" iJtE - E " t, "" , " +A)a 

| Wl _ W2 |-i/2^ 5 (5) 



hence 



When the temperature increases, the bounded vortices appear in increasing 
number. Their influence becomes more prominent, producing a deviation from the 
linear spin- wave contribution in equation (^), but the order parameter correlation 
function still decays algebraically with an exponent T]a-(T) which depends on the 
temperature. The transition eventually takes place at a temperature Tkt (usually 
called Berezinskii-Kosterlitz-Thouless critical temperature) when the pairs dissociate, 
thereby destroying short-range order in the system. The mechanism of unbinding 
of vortices was studied by Kosterlitz and Thouless jyj [l2|, [TsJ using approximate 
renormalization group methods. Above Tkt (in the vortex phase), the correlation 
function recovers a usual exponential decay. This very peculiar topological transition 
is characterised by essential singularities when approaching the critical point from the 
high temperature phase, t = T — Tkt ~~ * + , 

e~e^", (7) 

X~ (8) 

whence i) a = r\ a (Tkt ) = 2 - b x /b^. For reviews, see e.g. Refs. {J |||, [0| p[ 0. 

The existence of a scale-invariant power-law decay of the correlation function 
implies that at each value of the temperature there corresponds a fixed point, or 
equivalently that there is a continuous line of fixed points at low temperatures. The 
temperature is thus a marginal field in the spin-wave phase which becomes relevant 
in the vortex phase. 

Not much is known in the intermediate regime between the spin wave 
approximation at low temperature and the Kosterlitz-Thouless results at the 
topological transition and in fact the precise determination of the critical behaviour of 
the two-dimensional XY model in the critical phase remains a challenging problem. 
Most of the studies were dedicated to the determination of critical properties 
at Tkt t- Many results were obtained using Monte Carlo simulations (see e.g. 

| When approaching the transition from above, the value of Tkt is often obtained by fixing a = 1/2 
in equation (7), and ??ct(Tkt) then follows from a fit of the susceptibility to equation (8). 
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Refs. || n g | 
high-temperature series expansions 
existence of logarithmic corrections, e.g. 



m imp 

ansions |30|, |3l| , 1 



28, 29 ) at Tkt and slightly above, or 
but the analysis was made difficult by the 



X~£ 2 -*-(ln0 



211 



1 + 



In ln£ 
InT 



in the high-temperature regime, or 

(1d>i - w 2 |) 2e 



1 + 



In In |u>i — W2I 
In \wi - W2I 



with a probable value of 9 = -j^ 



(9) 



(10) 



Due to these 



at T KT exactly g g 
logarithmic corrections which make the fits quite difficult to achieve safely, the values 
of Tkt and 77,7 (Tkt) were a bit controversial as shown in table 1 of reference [p9[| . The 
resort to large-scale simulations was then needed in order to confirm this picture |28| . 
Although essentially dedicated to the KT point, some of these papers also report the 
value of r\ G (T) computed from the correlation function decay at several temperatures 
in the critical phase |L8| ^tj, ^2[ ^3) , while a systematic study of the helicity modulus 
at low temperature previously led more directly to the same information ]3E| . 

Recently, we proposed a rather different approach ^ which only requires a 
moderate computation effort. Assuming that the low-temperature phase exhibits all 
the characteristics of a critical phase with conformal invariance pq] (invariance under 
rotation, translation and scale transformations, short-range interactions, isotropic 
scaling), we use the covariance law of n— point correlation functions under the mapping 
of a two-dimensional system confined inside a square onto the infinite or half-infinite 
plane. The scaling dimensions are then obtained through a simple power-law fit where 
the shape effects are encoded in the conformal mapping, via the definition of a rescaled 
distance variable adapted to the description of the confined geometry. 

This technique is well known since the appearance of conformal invariance. 
Conformal mappings have been extensively used, mainly in the case of Ising or 
Potts models, in order to investigate the critical properties in restricted geometries. 
Density profiles or correlations have been investigated in various confined systems 
(surfaces p | B(f , c orners Q stripsjp], [l|, |4(J, squares Q Q or parabolic 
shapes f49|, 50[ |5l| for a review, see |j52fl). The exact expression of the two-point 
correlation functions have also been calculated on the torus for the pure Ising case |53| , 
and large scale Monte Carlo simulations of the 2d Ising model have been shown to 
reproduce the correct behaviour, including corrections due to the finite size of the 
lattice fH, ||. 

In this paper, we report a systematic application of conformal mappings to the 
investigation of critical properties of the 2c? XY model in the spin-wave phase. The 
physical quantities are obtained using Monte Carlo simulations in square shaped 
confined systems. We may then determine accurately the correlation function 
exponent rj a {T) in the low-temperature phase T < Tkt @ but also, by choosing 
convenient boundary conditions (BC), the surface critical exponent which describes 
the decay of the correlation function parallel to a free surface, ?7||(T) [[37J. The 
thermal exponent % will also be considered. Section|| presents the expected functional 
expressions of correlation functions and density profiles with various BC, and use of 
these expressions is systematically reported in section ^. It is worth noticing already 
that the results following the analysis of the two-point correlation functions are not 
fully convincing, but the main conclusions are then supported by the most refined 
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investigations of the density profiles. A discussion of the results is given in the last 
section. 

2. Functional expressions of correlation functions and density profiles in 
a square-shaped critical system 

2.1. Preliminary and notations 

In the following, we shall specify two different types of geometries, depending on our 
purpose. They are related through a conformal transformation. The critical properties 
in the first system will be reached through the study of the confined system in the 
second one. 

i) The infinite plane 

z = x + iy, 

— oo < Re z < +oo, 

— oo < Imz < +oo (11) 

or the upper half-plane 

z = x + iy, 

— oo < Re z < +oo, 

y = lmz>0. (12) 

This is the system which is eventually interesting in the thermodynamic limit. 
It is of course not directly accessible from numerical techniques. In the second 
case, the existence of a free surface materialised here by the real axis will allow 
to define both bulk and surface quantities. 

ii) The square 

w = u + iv, 

-L/2 < Rcw < L/2, 

0<lmw<L. (13) 

This is the natural geometry for Monte Carlo simulations. The boundary 
conditions may be open at the edges (square) or periodic (torus). In the following, 
we will obtain numerical results in the square or the torus, and our aim will be 
to translate them into their half- infinite or infinite system counterpart §. 

It is worth referring also to the infinitely long strip of width L, k + il, — oo < k < 
+oo, < I < L: this geometry is incidentally mentioned here because the choice of the 
boundary conditions is more obvious there than in the square system. Although no 
computation will be performed in this geometry, we will refer to it when introducing 
the profiles. It is particularly interesting, since many results concern such strips, e.g. 
when obtained with transfer matrices |5q ] or quantum chains [ |57| which correspond 

§ Here and in the following, z will constantly refer to a complex coordinate in the plane, while w will 
refer to the complex coordinate in the square geometry. In the plane, uhp specifies the upper half- 
plane. In the square, pbc, fbc, Fbc, and Ffbc means that the boundary conditions are respectively 
periodic, free, fixed, or mixed fixed-free according to a description which will be given later. 
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to the Hamiltonian limit of classical systems confined in strips. In the transverse 
direction (I— coordinate), the boundary conditions may be open (strip) or periodic 
(cylinder). A very accurate numerical study of the 0(n) model (including XY) using 
the transfer matrix of a related loop gas model was for example reported by Blote 
and Nienhuis ]5l| . The periodic strip is obtained from the plane through the standard 
logarithmic mapping k + il = In z which implies an asymptotic exponential decay 
for the correlation functions along the strip at criticality, 

(er kl+i ier k2+ ii) phc ~ exp(-27ra; cr |A;i - k 2 \/L). (14) 

The scaling dimensions are thus determined by the universal correlation length 
amplitudes, ^ 1 — 2irx a /L, which are given, in the transfer matrix formalism, by the 
two leading eigenvalues of the transfer matrix ^ 1 = ln(Ao/Ai). Blote and Nienhuis 
obtained the most precise values of the thermal and magnetic scaling dimensions at 
the KT transition, x e = \r] e = 2.000 000 (2) and x a = = 0.125 000 (1). In the 
following, we will exploit expressions similar to equation (|14|), transposed in the square 
geometry. We mention here that our approach proceeds from Monte Carlo simulations 
and cannot pretend to an accuracy comparable to the results of Blote and Nienhuis. 

2.2. Conformal rescaling of the correlation functions 

In an infinite system at criticality, the correlation function is known to behave 
asymptotically as 

(<x Zl .<x 22 )~ \z 1 -z 2 \- r >". (15) 

This defines the bulk correlation function exponent r\ a . The surface critical properties 
are also interesting and extra universal critical exponents can be measured from 
simulations performed in finite systems with open boundary conditions. In the semi- 
infinite geometry z — x+iy (the free surface being defined by the x axis), the two-point 
correlation function is fixed up to an unknown scaling function. Fixing one point z\ 
close to the free surface [z\ = i) of the upper half-plane (uhp), and leaving the second 
point Z2 explore the rest of the geometry, the following behaviour is expected: 

(cr Zl ■ cr Z2 } nhp ~ {yxy 2 )~ Xa i>{w), (16) 

with rja = 2x a . The dependence on 

ym 



zi - z 2 I 2 



(17) 



of the universal scaling function -0 is constrained by the special conformal 
transformation p9| , p0[ , and its asymptotic behaviour is implied by scaling e.g. 
ip(uj) ~ uj x " when t/2 > 1, with x\ = hrjn the magnetic surface scaling dimension. 

In a square system, the order parameter correlation function (cr Wl -cr W2 ) is strongly 
affected by boundary and shape effects and thus deviates significantly from the 
thermodynamic limit expressions (15) and (|T6|). It should obey a scaling form which 



reproduces the expected power-law behaviours, for example 

(<r Wl ■ cr W2 ) = \wi-w 2 \~ TI "f s < l Xwi/L,W2/L), (18) 

where the function / sq . encodes shape effects in a very complicated manner which also 
depends on the boundary conditions, so that the connection between equations (|15l) or 
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( |l6| ) and (|l8|) is not obvious. Conformal invariance provides an efficient technique to 
avoid these shape effects, or at feast, enabfe to inciude explicitly the shape dependence 
in the functional expression of the correlators through the conformal covariance 
transformation under a mapping w(z): 



(19) 



Conformal invariance holds in isotropic systems with short range interactions which 
exhibit translation, rotation, as well as scale invariance. Usually these conditions are 
so restrictive that they can only be valid at the critical point, and we assume here 
that the general covariance equation ( |l9| ) holds in the whole T < Tkt phase. 

In order to get a functional expression of the correlation function inside the square 
geometry w, one simply has to use the convenient mapping w(z): 

i) Conformal mapping of the plane onto the torus: in the continuum limit, the 
correlation function of the Ising model on the torus is exactly known |5^, [h], [58| 



r W2 ) pbc 7 



Ising 



E 10,(0)1 



y=2 



eK( 



Wl2 

2L 



1/4 



(20) 



(rj a — 1/4) where 6 v (a) are the Jacobi theta functions and w\2 — W\ — u>2- 
A similar expression exists for the energy-energy correlations fllPf . It can be 
generalised to a critical system characterised by any correlation function exponent 
r\ a . For simplicity, we fix in the following one point (wi) as the origin, and the 
second point is simply written w. The dependence on the relative distance w in 
the above expression, 



El^ (w/2L)\ x \9 X (w/LT 1 ^, 



(21) 



is compatible with the general covariance equation (|l9|), rewritten as 

(to • v w ) phc - \w'(z)\-^/ 2 [z(w)}-^, (22) 

provided that |^(0)| -C \z(w)\ (large distance behaviour). The r.h.s. is a function 
of w hereafter denoted by [C{w)]~ v " , with r] a — 1/4 and the rescaled distance is 



aw) = \9 1 (w/L)\x El0»K2£)l 



(23) 



Therefore equation ( |20| ) can be generalised to any correlation function of a critical 
system inside a torus, 



r w) pbc 



(24) 



A sketch of the correlation functions in the plane and in the square geometries is 
shown in figure [l]. 



Figure 1. Correlation functions in the plane (left) and in the square (pbc) geometries 
(right). 



) Conformal rescaling of the upper half-plane inside a square: the conformal 
transformation of the half-plane z — x + iy (0 < y < oo) inside a square w — u + iv 
of size LxL(—L/2<u<L/2,0<v<L) with open boundary conditions along 
the four edges is realized by a Schwarz-Christoffel transformation Kfl 



2Kw 



(25) 



(26) 



sn(2Kw/L) the Jacobian elliptic sine, K = K(k) = F(l,k) the complete elliptic 
integral of the first kind, and the modulus k depends on the aspect ratio of A and 



Here, F(z, k) is the elliptic integral of the first kind, 
F(z,k)= / [(l-i 2 )(l-fc 2 t 2 )]- 1 / 2 dt, 



is here solution of K(k)/K(\/1 — k 2 ) 



k = 4 



0.171573, q = c~ 



(27) 



Using the mapping (p5|), one obtains the local rescaling factor in equation (|l9|), 
w'(z) = ^[(1 — 2 2 )(1 — fc 2 z 2 )] -1 / 2 , and inside the square, keeping wi ~ 0(1) 
fixed, the two-point correlation function becomes (see e.g. Ref. »□)) 



k(w) = Im 



2Kw~ 



L 



x 



1 - SI 



,2KuA 



I - k 2 sn 2 



,2Kw 



-1/2 



(28) 



where fbc specifies that the open square has free boundary conditions. This 
expression is correct up to a constant amplitude determined by k(wi) which is 
kept fixed, but the function ip(u>) is still varying with the location of the second 
point, w. 

A sketch of the correlation functions in the plane and in the square geometries is 
shown in figure ^| 



Figure 2. Correlation functions in the upper half-plane (left) and in the square (fbc) 
geometries (right). 



2.3. Density profiles in restricted geometries 

In order to cancel the role of the unknown scaling function in the half-plane geometry 
in equation (28), it is more convenient to work with a density profile m(w) in the 
presence of symmetry breaking surface fields h^A on the boundary <9A of the lattice 
A. This is a one-point correlator which scales in the half-infinite geometry as: 

m(z) = (a z ■ h aA(2) ) uhp = constx?/" 2 ^ (29) 

and it maps onto (see figure [s]) 

m Fhc (w) = (er w • h aA ( w) ) F bc = const x[k(w)]~% v ° (30) 

where the function k(w) defined in equation (28) again comes from the mapping and 
Fbc means that the open square has fixed boundary conditions (figure ||) . 

Equation ( p9| ) is a special case of conformally invariant profiles predicted e.g. by 
Burkhardt and Xue [H, [45| in strip geometries k+il or in their half-infinite counterpart 
z = x + iy: 

m ab (z)=y- x °F ab (cos6). (31) 

The scaling function F ab (x) depends on the universality class and on the boundary 
conditions denoted by a and b pM . This notation is reminiscent from the strip 
geometry where a and b specifies the two boundaries of the strip at I = and I = L 
respectively. In the case of fixed BC (denoted by Fbc in the square and FF in the 
strip) corresponding to equation ( p9| ) and (|3^) , one gets 

m Fhc (z) = <T W2 , (32) 

while in the case of fixed-free BC (Ffbc in the square or more naturally Ff in the strip) 
|| Burkhardt and Xue have shown in the case of Ising and Potts models that 

m Ffbc (z)-^/ 2 (cosi0)""/ 2 . (33) 

|| The spins located on the x > half axis (1 = strip surface) are kept fixed while those of the x < 
half axis (I = L strip surface) are free. 
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Figure 3. Order parameter profile in the plane (left) and in the square geometries 
(right) with fixed boundary conditions (top) and mixed fixed- free BC (bottom). 



This expression was since then mainly used in the strip geometry |32|, p3fl . In the 
square geometry, the fixed- free BC correspond to keeping the spins fixed for u > 0, 
and the profile is expected to obey the following ansatz |57| (figure ||) : 



TOFfbc(w) ~ const x [k(w)] [n(w)} 2,711 , 



M«0 = ^ [ 1 + Rc 



2Kw 



L 



2Kw 



In this paper we will essentially apply equations (g4|), (28), 
to determine r\ a and 7711 . 



(34) 

and (34) in order 



3. Monte Carlo simulations 

3.1. Description of the algorithm 

Simulations of 2d XF-spins are performed using Wolff's cluster Monte Carlo 
algorithm |64|. The boundary conditions may be periodic, free or fixed. When 
the boundaries are left periodic or free, there is no particular problem, but a 
relatively precise numerical determination of the order parameter profile of the 2d 
XY model confined inside a square with fixed boundary conditions (playing the role 
of ordering surface fields hgA(«;)) introduces a priori a technical difficulty which is 
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easily circumvented. In practice, the symmetry is broken by keeping the boundary 
spins located along the four edges of the square fixed, e.g. er w — (1,0), Vu; G <9A 
during the Monte Carlo simulation. The Wolff algorithm should thus become less 
efficient, since close to criticality the unique cluster will often reach the boundary 
and no update would be made in this case. To prevent this, we use the symmetry of 
the Hamiltonian (2) under a global rotation of all the spins. Even when the cluster 
reaches the fixed boundaries dA(w), it is updated, and the order parameter profile 
is then measured with respect to the common new direction of the boundary spins, 
rriFbc(w) — (cr w ■ <TaA(u)))Fbc- The new configuration reached would thus correspond - 
after a global rotation of all the spins of the system to re-align the boundary spins in 
their original (1, 0) direction - to a a new configuration of equal total energy and thus 
the same statistical weight as the one actually produced. 

In the literature, many papers were dedicated to the study of the characteristic 
properties of the Kosterlitz-Thouless transition. The value of Tkt will be considered 
as known with a good accuracy %. We will take here the value reported by Gupta 
and Baillie ||, or Campostrini et al @, k B T KT /J = 0.893(1). 

3.2. Behaviour of the correlation functions 

3.2.1. Scaling of correlations on the torus: Simulations are performed on a square 
with periodic boundary conditions. Here we choose systems of size L = 100. After 
thermalization (10 6 iterations) and computation (10 6 other iterations), we get the spin 
orientations {0 W } at all the lattice sites. The correlation function is defined as 

(er • cr w ) phc = (cos(0 o - w )) pbc , (35) 

where is an arbitraty reference site (e.g. a corner in figure Q). 

It has to be fitted to equation ( p4| ) with C( w ) given by equation (|2^). Expansions 
of the 0- functions can be found in the literature, e.g. in Abramowitz and Stegun |65fl, 



O^w/L) = 2e-' r / 4 ^(-l) n e- n (™ +1 ) 7r sin 

71=0 

and similar expression for 02, 03 and 0\ . Here it has been taken into account the fact 
that the system is inside a square (not a rectangle L x Lf), therefore the expansion 
parameter e~ 7rL l L reduces to q — e _7r . Keeping only the leading terms, we obtain 

0i(w/L) = 2e- 1T/4 (sm7rw/L-e- 27T sin Zirw/L + ...) (37) 



(2n + l)nw 



(36) 



4 

\0„(w/2L)\ = 2(1 + e-* /4: (smiTw/2L + costtw/2L) 

- e- 97r/4 (sin37ru;/2L - CQ^wj2V) + ...). (38) 
Using w = u + iv and defining the symbols 
sc[u,w] = sin(7ra/-L) cosh(nv/L) 

cs[u, v] = cos(iru/L) s'mh(irv/L) (39) 

If The numerical value of the transition temperature has been deduced from different approaches, e.g. 
high-temperature series expansions (e.g. in El), temperature dependence of the correlation length 
(e.g . in ]23|), zeroes of the partition function (e.g. in J29|), jump in the helicity modulus (e.g. in 
|20j). . . . and the recent estimations are compatible with each other. 
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Figure 4. Log-log plot of the order parameter correlation function in a square 
system with periodic boundary conditions vs the rescaled variable ((w). The data 
correspond to a quarter of the lattice A, for < u < L/2 and < v < L/2, L = 100. 
Four different temperatures in the spin-wave phase are shown (kgT/J = 0.1, 0.3, 
0.5, and 0.7 from top to bottom), and at the KT transition temperature. 



and similar expressions for ss[u,w] and cc[u, v] which enable a compact notation for 
shnrw/L — sc[u, v] +i cs[u, v], the expression of C( w ) m equation ( p3| ) is given in terms 
of trigonometric functions: 

\9i{w/L)\ = 2q [(sc[u, v] - q 8 sc[3w, 3v]) 2 + (cs[u, v] - q 8 cs[3w, 3v]) 2 ] 1/2 

(j2 \0uH2L)^j = 2 [1 + g( S c[|, f ] + cc[|, f ]) - g 9 (sc[f , f ] - cc[f , f ])] 2 

+2 [1 + g(cs[f , §] - ss[|, §]) - g 9 (cs[f , f ] + ss[f , f ])] 2 (40) 

It is instructive to recover the cylinder limit from the torus geometry. For a square 
of arbitrary aspect ratio L' /L, the expansion parameter q in the previous series 
expansions has to be replaced by e _7ri / L , and the larger the aspect ratio, the faster 
the series converges. Let us consider the limit L' — * oo, L fixed. Keeping only the 
leading terms and omitting unimportant prefactors, we have 

\0i(w/L)\ - \sc 2 [u,v}+cs 2 [u,v}\ 1/2 , (41) 

TJ 1 1 'j— »nn 
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hence from C( w ) ~ \6i( w /L)\ it follows the traditional expression |36| 

(42) 

which leads to equation (|lj) in the limit v/L ^s> 1 + . 

We can then plot the correlation function data directly versus the rescaled variable 
£(u>) on a log-log scale, as shown in figure ||. On this scale, a linear behaviour would 
indicate a power-law decay. A deviation from the straight line can be observed 
at large values of ((w) where the data are scattered. The quality of the data is 
unfortunately not improved by a better statistics (ten times more iterations). This 
correction is more pronounced as the KT point is approached where one knows that 
logarithmic corrections are present. From this figure, it is clearly difficult to get a 
reliable estimation of the exponents, but we nevertheless performed a fit in the range 
£(u>) < 4.10 -4 where the behaviour is closer to a pure power law *. The straight lines 
reported there correspond to (power-law fits with) slopes ri a {T) as given in table [l] 
in the conclusion. They match perfectly the correlation function data in the regime 
of small (,{w), therefore if they do not confirm nicely the expression i n (p4h , they at 
least do not contradict this equation. A stronger support to equation fl24|) is the fact 
that the values of the exponents reported in table [l] are consistent with more refined 
estimates deduced from the fit of the order parameter profile in section |3.3| (the reader 
can compare columns a and b in table [I] to convince himself of the compatibility of 
the results). 



{(TO ■ (T m/pbc 



cosh (_j_ cos /_ 



3.2.2. Scaling of correlations in open square: Simulations are then performed 
according to a similar procedure, but now in an open square with free BC along 
the four edges. Equation (28) may be used to estimate the surface exponent rj\\(T). 
Plotting the quantity 



/IT). 



x [«(«;)] ^'CH =vh 



as a function of the variable u defined in equation (|1 
w-coordinates as follows: 



(43) 

and given in terms of the 
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sn — ; sn 



L 



(44) 



will give access to the universal function ip(cu). Here, the values of r] a (T) are taken 
from the previous section (column a in table |l|). As we already mentioned above, 
these are not our most precise determinations of the bulk scaling dimensions, but 
the deviation with the forthcoming results of the following section being small (see 
table |l|) , there is no dramatic consequence to use them. 
The function ip(oj) is expected to behave as 



(45) 



+ In the standard notation, the roles of u and v are reversed, because we have considered here the 
infinite direction of the cylinder in the imaginary direction. 

* This range does not correspond to the right limit for an investigation of the asymptotic behaviour 
of the correlation function (especially when logarithmic corrections may be involved) , since the limit 
\z\ — » oo in the plane corresponds to C( w ) ~* °°- 
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Figure 5. Log-log plot of the universal scaling function ip(u) at the KT point and 
below. The dashed lines represent a fit leading to the surface correlation function 
exponents (L = 101). 



and leads to a determination of the surface exponent as a function of the temperature. 
In the case of the 2d Ising model for example, the scaling function was calculated 
analytically by Cardy |39|, 



Rising M - ^(l+4^) 1 /4 + (l + 4^-1/4 



// 2 (1-C+ V +...), 



(46) 



in accordance with Jj|, sms = 1. 

The rescaled space variables k(w) and u> involve Jacobi elliptic sine. The series 
expansion is the following, 



2Kw 2tt 



sn ■ 



\ " 



-7r(n+l/2) 



L 



kK ^ 1 - e-^n+i] 

n— 



sin 



(2n + l)nw 



L 



(47) 



and the plot of ^{oj) is shown in figure at several temperatures below and at the 
KT transition. The surface scaling dimension ^n\\ (T) is deduced from a power-law fit 
at small values of u> < 0.03 (dashed lines) j. At the Kosterlitz-Thouless transition, we 



tj, This is the right limit this time, since | jz | — > co corresponds to u> — > 0. 
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can also fit the scaling function for comparison with the Ising case. One obtains 

V-ktM ~ w a274 (l - 0.401cj + 0.221w 2 + ...). (48) 

3.3. Behaviour of the density profiles 

3.3.1. Scaling ofmFbc(w): Density profiles present the advantage of being one-point 
functions. They are thus supposed to be determined more precisely numerically than 
two-point correlation functions and should therefore lead to more refined estimations 
of the critical exponents. We start with simulations inside an open square with fixed 
boundary conditions (Fbc), according to the prescriptions given in the description 
of the algorithm. The order parameter profile mpbc(w) is defined according to the 
definition of a reference orientation on the fixed boundaries, i.e. 

mFbc{w) = (cr w ■ CTdA( w )}Fbc 

= {COS(0 W - 03A(xt>)))Fbc- (49) 

After averaging over the 'production sweeps', one gets a characteristic smooth profile 
as shown in figure ^ which also presents a typical configuration with fixed BC. 




Figure 6. Monte Carlo simulations of the 2d XY model inside a square 48 X 48 spins 
with fixed boundary conditions below the Kosterlitz-Thouless transition temperature. 
Typical configuration on the left (where the fixed BC are easily observed by the 
common direction taken by all the border spins) and average over 10 6 MCS/spin 
after cancellation of 10 6 for thermalization (cluster update algorithm) on the right. 



A log- log plot of mpbc(w) with respect to the reduced variable k(w) is shown in 
figure ^ at two different temperatures below Tkt, roughly at Tkt, and one temperature 
slightly above. A first observation is the confirmation of the functional form of 
equation ( |30|) in the whole low-temperature phase. One indeed observes a very good 
data collapse of the L 2 points onto a single power-law master curve (a straight line on 
this scale). The next information is the rough confirmation of the value of the critical 
temperature, since above Tkt, the master curve is no longer a straight line, indicating 
that the corresponding decay in the half-infinite geometry differs from a power-law 



as it should in the high-temperature phase ff. In reference 36 , we have shown how 



ff One should nevertheless mention that this does not lead to a precise determination of Tkt, since 
at kgT/J = 0.9 for example, the master curve is hardly distinguishable from a straight line. 
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Figure 7. Monte Carlo simulations of the 2d XY model inside a square 100 X 100 
spins (10 s MCS/spin after cancellation of 10 6 for thermalization, cluster update 
algorithm). The figure shows the rescaled order parameter density against the 
rescaled distance on a log-log scale at different temperatures below the Kostcrlitz- 
Thouless transition temperature and slightly above. 



the computation the x 2 P er d.o.f. as a function of temperature is an indicator of 
the location of the KT transition. The \ 2 has a very small value, revealing the high 
quality of the fit, in the low-temperature phase and then increases significantly above 
Tkt when the behaviour of the density profile is no longer algebraic. 

We have also shown that there is probably no logarithmic correction at the 
Kosterlitz-Thouless point, since the curve of figure |7| does not seem to differ 
significantly from a true power-law. Although not a proof, a numerical evidence of 
that is the plot of WFbc(w) x [k(ui)] 1 / 8 vs k(w) which remains constant on the whole 
range of values of k(w) [ p6| . This observation makes the order parameter profile with 
fixed boundary conditions a very convenient quantity for the determination of the 
bulk correlation functions exponent, since it displays a pure algebraic decay up to the 
Kosterlitz-Thouless transition point f . 

Fitting the curves of figure ^ to power-laws leads to the values of the scaling 
dimension ri a (T) in the whole critical region T < Tkt as shown in figure 10 in the 
conclusion. The data are also reported in table |l| in the conclusion. The accuracy is 



f The decay exponent, if not fixed, but determined numerically, leads to a quite good result (for 
example for a size L = 100, we get rj a (Tkt ) = 0.250(28)). 
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much better than in the case of the two-point correlation functions. 

3.3.2. Scaling of mFfbc(w): For this last choice of boundary conditions, hereafter 
referred to as "fixed- free BC" (and denoted by Ffbc), the open square has half of its 
boundaries which are left free while the other half is kept fixed. This is illustrated in 
figure ^ where the averaged profile is also shown at two different temperatures. The 
shape of the plotted surface helps understanding the BC. 




Figure 8. Monte Carlo simulations of the 2d XY model inside a square 48 X 48 spins 
with mixed fixed-free boundary conditions below the Kosterlitz-Thouless transition 
temperature. Typical configuration on the left and average over 10 6 MCS/spin after 
cancellation of 10 s for thermalization (cluster update algorithm) at ksT/J = 0.488 
and k B T/J = 0.915 on the right. 

In the previous section, we reported some evidence that the logarithmic corrections 
at the KT point are negligible for the order parameter profile with fixed BC. Assuming 
that there will be no substantial difference for fixed-free BC, we perform a fit of 
lnmFfbc(w) against the two-dimensional linear expression, const — ^r] a (T) In. k(w) + 
5?7||(T)ln/i(it;), as expected from equation (34). The value obtained for the bulk 
correlation function exponent ^(T) (already known precisely from previous section) 
will be a test for the reliability of the fit. The simulations are performed at several 
temperatures, up to ksT^T/ J ~ 0.893. The three parameters linear fits lead to correct 
determinations of both exponents (we studied different system sizes ranging between 
L = 48 and L — 100 and even 200 at the KT point). The value of the bulk exponent 
is in a good agreement with the determination which follows from the fit of the profile 
with fixed BC, and the surface exponent also fits nicely to the values deduced from 
the two-point correlation function with free BC. 
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3.3.3. Scaling of Aepbc-fbcf^)-' We concentrated mainly up to now on the magnetic 
properties. The thermal exponents are also interesting quantities which characterise 
a critical point. In the case of the XY model, the temperature is a marginal field, 
responsible for the existence of a critical line in the whole low-temperature phase. It 
thus implies a thermal scaling exponent x £ — d — y t = 2 which ensures a vanishing 
RG eigenvalue y t = (up to Tkt where it is consistent with the essential singularity 
of £ above the KT point). The energy-energy correlation function should thus decay 
algebraically as 

(e Zl s Z2 ) ~|«i-z a |-"«, (50) 

with r) s (T) = 2x £ = 4 VT < T KT . 

We check qualitatively this expression of the energy-energy correlations from the 
behaviour of the energy density profile. The energy density at site w is for example 
defined as the average value of the energies of the four links reaching w. 

£w = 7T7 y~]( <T «'-/i ' & w + & w ' &w+p)- (51) 




Figure 9. Monte Carlo simulations of the 2d XY model inside a square 100 X 100 
spins (10 s MCS/spin after cancellation of 10 7 for thermalization, cluster update 
algorithm). The main figure shows the local energy density vs the rescaled variable 
k(w) at different temperatures and for different BC. The insert is a log-log plot of 
the difference Ae(ui) as explained in the text. 

In comparison with the magnetisation profile, a difficulty occurs due to the 
existence of a regular contribution in the energy density. This term cancels after 
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a suitable difference between profiles obtained with different BC. For that reason, we 
made two different series of simulations with free (fbc) and fixed (Fbc) boundary 
conditions. From the latter simulations we had already extracted the singular 
behaviour of the magnetisation profile since there is no regular contribution in the 
scale-invariant phase, 

m Fbc (z)=A m (T)y~^ T V 2 , (52) 

where A m (T) is a non-universal "critical" amplitude f which depends on the 
temperature. 

On the other hand, the two series of simulations are necessary in order to extract 
the singularity associated to the energy density which contains a non universal regular 
contribution (so(T)} which depends on T and a singular contribution: 

(e,) F bc = (eo(T)) + B Fhc (T)y-^ [T)/2 , (53) 
(e,)fbc = (eo(T)) + 6 fbc (T)y-"-( T )/ 2 . (54) 

This is clearly illustrated in figure ^ where convergence towards the same temperature- 
dependent constant (eo(T)) is shown. We also observe that the amplitudes of the 
singular terms have opposite signs. Therefore a simple difference of the quantities 
measured in the square geometry, 

Ae(w) = (e w ) Fhc - (e^fbe ~ AB x [ K (w)]-i^ (55) 

leads to the value of the thermal scaling dimension r] s (T). 

The insert in figure ^ shows a log-log plot of the difference Ae(w) vs k(w) at three 
temperatures below Tkt- Due to strong fluctuations, the data scatter around straight 
lines which represent the theoretical slopes [k(w)]~ 2 . This figure, though far from 
being definitely conclusive, confirms that the exponent of the decay of energy-energy 
correlations keeps a constant value r\ e (T) = 4 in the low-temperature phase of the XY 
model. At the KT transition, this value was confirmed very accurately by Blote and 
Nienhuis f36|] . 

4. Discussion 

The different determinations of the bulk and surface correlation function scaling 
dimensions reported in this paper are mutually consistent. These exponents are 
obtained using completely independent techniques, from four independent series of 
simulations (the boundary conditions, but also the quantities which are computed 
are different). The fits themselves are different, since these quantities obey specific 
functional expressions. 

The exponents are plotted in figure [l(] for different system sizes. Our results (for 
a few values of the temperature) are also given in table |l|. In column a, 6, and c we 
present the bulk exponent. Columns d and e give the surface exponent. The error 
bars (for the values resulting from the fits of the magnetisation profiles, columns b, c, 
d) correspond to one standard deviation. We have to mention here that these error 
bars should be considered with some care for at least two reasons. A first reason 
is that we did not make any extrapolation to the thermodynamic limit, because the 
results at a few sizes seem already very stable as it can be observed in figure |io| . 

\ Here we use the term "critical" , since the whole low temperature phase is considered as critical. 
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Figure 10. Bulk and surface correlation decay exponents of the 2d XY model as 
a function of the temperature. The values follow from independent linear fits of the 
magnetisation profiles with Fbc (circles) or Ffbc (triangles). 



Although we did not analyse systematically the size effects, we believe that this is 
not the most important contribution to the error. A second reason has more drastic 
consequences. This is the influence of the fitting window. When the range of fit is 
varied in a reasonable scale in the 'large k' regime, it usually produces a small deviation 
in the resulting exponent. For example, fitting a power law ~ [^(w)]^ in the window 
0.10 < k(w) < 1 or 0.25 < n(w) < 1, lead to compatible results within the error 
bars (in this example at the KT transition, we get respectively r\ a = 0.260(18) and 
0.250(28)). If the window range is changed significantly on the other hand, we do no 
longer control the stability of the results, since the errors increase considerably (for 
example r\ a = 0.247(92) for 0.50 < k(w) < 1). Because of these constraints, we believe 
that the values of the exponents that we measure are reliable, but the errors are only 
rough estimates. 

The weaker result of the paper probably comes from the correlation function 
with periodic boundary conditions which was fitted to an approximate functional 
expression. There are also in this case logarithmic corrections and an unavoidable 
scattering of the data which prevent a fit in the interesting asymptotic region. The 
consequence is a bulk exponent which is a priori obtained with a poor level of 
confidence, and for this reason we decided to avoid any error bar in column a of table [j]. 
It is amazing to observe that in spite of these restrictions, these results do indeed agree 
with the most refined values deduced from the magnetisation density profile with fixed 
boundary conditions (see column b in table |l|) . Another consequence of our approach 
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Table 1 . Values of the exponents of the bulk and surface correlation function in the 
spin-wave phase of the 2d XY model. They are obtained after fitting the density 
profiles or correlation functions to the conformal functional expressions for a system 
of size L = 100. We have checked that the results are already very stable at such a 
size. We mention that determinations a, b and c are completely independent for the 
bulk exponent, and that determinations d and e for the surface exponent are also 
independent. For determination e, the value of the bulk exponent in a is used as 
input. The most reliable results are printed in bold font. 



Correlation function exponents 








b 


c 


sw. 


d 


e 


SW. 


0.1 


0.017 


0.017(2) 


0.017(1) 


0.016 


0.032(1) 


0.031 


0.032 


0.2 


0.037 


0.035(4) 


0.036(3) 


0.032 


0.069(1) 


0.064 


0.064 


0.3 


0.057 


0.052(7) 


0.052(5) 


0.048 


0.103(2) 


0.100 


0.095 


0.4 


0.077 


0.073(9) 


0.074(6) 


0.064 


0.143(4) 


0.140 


0.127 


0.5 


0.103 


0.096(11) 


0.100(8) 




0.186(5) 


0.181 




0.6 


0.127 


0.122(14) 


0.122(13) 




0.231(8) 


0.228 




0.7 


0.154 


0.153(18) 


0.155(14) 




0.298(12) 


0.294 




0.8 


0.194 


0.188(22) 


0.187(29) 




0.374(17) 


0.366 




0.893 


0.249 


0.250(28) 


0.250(34) 




0.542(22) 


0.548 





a deduced from the fits of (ctq ■ °"w) P bc (L = 100), section 3.2.1. 
° deduced from the fits of mpi, c (to), section 3.3.1. 

c deduced from the fits of cnpfhcC'") (L = 100 or 200 at the KT point), section 3.3.2. 
™ deduced from the fits of mpfbc(l") [L = 100 or 200 at the KT point), section 3.3.2. 
e deduced from the fits of (a wi ■ &w)fbc {L = 101), section 3.2.2. 
S.W. gives the result of the spin wave contribution. 

is again the a priori weak level of confidence of the surface exponent deduced from the 
correlations in the open system (column e in the table), since we use the bulk values 
from column a as input parameters for this latter fit. The surface exponent in fact 
fits correctly the most refined values deduced from the magnetisation density profile 
with mixed fixed- free boundary conditions (column d) . Obviously, the stronger results 
come from the magnetisation profiles analyses. We believe also that another result 
of interest is the confirmation that the functional expressions derived in section 2 are 
valid and are in fact very efficient for the numerical investigation of critical properties 
of two-dimensional systems. Such applications of conformal mappings are of course 
well known in the case of strip geometries with e.g. the celebrated correlation length 
amplitude-exponent relation ]G6|. 

Let us now discuss briefly the results. In the low-temperature limit, one knows 
that the spin-wave description captures the essential of the physics of the problem. 
From this limit, the linear behaviour of the bulk exponent with the temperature, 
rja(T) = /csT/27rJ, is easily recovered. It is of course reasonable to expect that the 
low-temperature limit of the surface exponent can also be recovered by the spin wave 
approximation (see appendix). Indeed, the surface properties of the Gaussian model 
at the ordinary transition were studied by Cardy p9| |. They describe to the low- 
temperature regime of the XY model. The surface-bulk correlation function exhibits 
an algebraic decay 

/ irr \ — i/irK 

<cos(# -#r))uhp^ (— ) (56) 
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and the corresponding surface exponent is thus given in the spin wave approximation 
by the linear behaviour 

V\\(T) = ^-, T^O. (57) 

This result is simply twice the bulk value, rj\\(T) = 2r) a (T). This relation fits the 
numerical data at low temperature as shown in figure [T^ or in the table. 

At the KT transition, the bulk exponent is fully coherent with the result of 
Kosterlitz and Thouless RG analysis ^(Tkt) = 1/4. For the surface, the value 
?7||(2kt) — 0.54 is close to the known result | j39j, a value which cannot be excluded 
from our analysis (though out of the error bar), since the uncertainties that we 
quote are only indications. It is also coherent with the surface exponent obtained 
at the KT point of the quantum Ashkin- Teller chain (e = — 2 -1 / 2 in reference fj?!). 
t]||(Tkt) = 27r~ 1 arccos (— e) = 1/2. This value is furthermore in agreement with the 
relation r]\\{T) — 2rj lT (T) 1 provided that it is valid up to the KT point. This relation 
was established in the spin-wave approach, and to extend it to the vortex dominated 
regime, closer to Tkti we note that the effect of vortices (when T — > ?kt) when treated 
in the Coulomb gas approach [l^] is only to increase the effective temperature 
in the spin wave approximation, hence to enhance the decay of the correlations. 
Like ?7||(r) = 2rj a (T) = ksT/nJ is correct in the spin wave treatment, the relation 
V\\{Teft) = 2?7 CT (r o ff) = ksTcg/nJ should hold in the whole critical phase, with a shift 
in temperature which depends on the vortices interactions, 

fcsTcff = k B T - ivr 2 r 2 {q(0)q(r)), (58) 
|r|>o 

where (q(O)q(r)) are the vortex intensity correlations (which are negative in the 
neutral low temperature phase). It is thus reasonable to expect that the relation 
r^ll (T) = 2t7 ct (T) is exact for the whole critical phase of the model. It can be checked 
that the numerical values in table [l] support this expression. 

The approach of the KT point, where the temperature dependence of the bulk 
correlation function exponent is known, is also interesting. From the numerical data, 
r]a{T) can be expanded close to the KT point T — ► T^ T , showing a leading square root 
behaviour JlJ, as illustrated in figure |ll|. The same behaviour is obtained for the 
surface exponent. One should nevertheless mention here that the plot of (r) a (T) — j) 2 
or (t7|| (T) — i) 2 against T seem to indicate a slightly lower value of Tkt than 0.893. A 
clarification of that point is beyond the scope of this paper, since it would definitely 
require simulations of really large systems in order to rule out the possibility of finite- 
size dependence of the measured quantities. 

Eventually, we mention that our analysis of the energy density profiles is 
compatible with the marginality condition of the temperature in the low temperature 
phase of the XY model. This condition, equivalent to the continuous variation of the 
magnetic scaling dimension, indeed implies that the thermal scaling dimension keeps a 
constant value x e = 2 which do not depend on the temperature, and this value enable 
to fit reasonably the numerical data. 

In this paper, we performed an accurate study of the low-temperature phase of 
the 2d XY model using efficient techniques, based on a conformal mappings of the 
correlation functions and density profiles. At the Kosterlitz-Thouless transition, the 
results are in agreement with the general picture of the RG analysis and the critical 
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A L = 100 Fflyc 




Figure 11. Square root behaviour of the bulk and surface correlation function 
exponents as the Kosterlitz-Thoulcss point is approached. 

exponents describing the algebaric decay of the correlations in the low temperature 
phase are measured. 
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Appendix 

In this appendix, we recall the surface-bulk correlation function obtained in the spin 
wave approximation. The correlation function between a spin close to the surface 
(ctq) and another spin far in the bulk (<x 2 ) is defined according to (<r • <x 2 ) u hp = 
(cos(#o— #z))uh P , where the subscript uhp means here that Imz > 0. We use a standard 
notation r for z from now on. As it is well known, in the harmonic approximation, the 
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Hamiltonian becomes quadratic and the thermal average leads to the Gaussian model 
which implies that 

{cos(9 - r )) uhp = e -5<( e o-0r) 2 >uh P _ (1) 

In Fourier space, 9 r = N~ x / 2 J^q exp(— iqr)#q, (N = L 2 ), the Hamiltonian becomes 
(the spins are located on the vertices of a square lattice of spacing unit a) § 



const — T-ft'^ 2(2 — cosaq x — cos aq y )\9q\ 2 



H 



k B T 4 



q 

const- Jif ^a 2 |q| 2 |0 q | 2 (2) 



q 

an leads to a partition function 

Hence, the quadratic fluctuations of the Fourier amplitudes are linear in T, 

<N 2 > = (4) 

and the integral (N' 1 £ q -» J dq) 

((^0-^r) 2 ) = |^ | dq(l-cosqr)(|0 q | 2 ), (5) 
is easily evaluated, leading to 

1 rl a dq f 27T 
((do - r ) 2 ) = -g— I —I d<p(l - cos qr cos if) 
I" K Jo 1 Jo 

nK J u 

In-, (6) 
nK a 

where a cutoff <~ it /a was introduced. The surface-bulk correlation function follows: 

/7rr\ -i/ 7r -ff 

(cos(0 o - 6 r )U P * (— ) (7) 

and the corresponding surface exponent is thus given in the spin wave approximation 
by the linear behaviour 

Vll(T) = ^f, T^O. (8) 

§ In contrast with the standard calculation in the whole plane, an extra factor i appears from the 
orthogonality relation in the uhp. 
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